A shift in lung macrophage composition is associated with COVID-19 severity and recovery

Although it has been more than 2 years since the start of the coronavirus disease 2019 (COVID-19) pandemic, COVID-19 continues to be a worldwide health crisis. Despite the development of preventive vaccines, therapies to treat COVID-19 and other inflammatory diseases remain a major unmet need in medicine. Our study sought to identify drivers of disease severity and mortality to develop tailored immunotherapy strategies to halt disease progression. We assembled the Mount Sinai COVID-19 Biobank, which was composed of almost 600 hospitalized patients followed longitudinally through the peak of the pandemic in 2020. Moderate disease and survival were associated with a stronger antigen presentation and effector T cell signature. In contrast, severe disease and death were associated with an altered antigen presentation signature, increased numbers of inflammatory immature myeloid cells, and extrafollicular activated B cells that have been previously associated with autoantibody formation. In severely ill patients with COVID-19, lung tissue–resident alveolar macrophages not only were drastically depleted but also had an altered antigen presentation signature, which coincided with an influx of inflammatory monocytes and monocyte-derived macrophages. In addition, we found that the size of the alveolar macrophage pool correlated with patient outcome and that alveolar macrophage numbers and functionality were restored to homeostasis in patients who recovered from COVID-19. These data suggest that local and systemic myeloid cell dysregulation are drivers of COVID-19 severity and modulation of alveolar macrophage numbers and activity in the lung may be a viable therapeutic strategy for the treatment of critical inflammatory lung diseases.

centrifugation at 300 rcf for 15 mins. Cell viability and counts were assessed by acridine orange and propidium iodide (AOPI) staining in automated Nexcelom Cellometer Cell Counters. PBMC were resuspended at a concentration of ~10x10 6 cells/mL in Human Serum and 10% dimethyl sulfoxide (DMSO) and stored at -80°C for 24 hours before transfer to liquid nitrogen storage.

Olink measurements of COVID-19 serum, data normalization, and clustering analysis
Olink was performed on COVIDvolunteer and COVID + patient serum samples in BSL2 + according to manufacturer instructions. Count (Ct) values were generated by Olink NPX manager software. To control for technical variability between plates, we included in each plate 2 technical control replicates from a single mixture of pooled blood from healthy donors and estimated a control value per plate, defined as: Where: a is a given analyte, i is a given sample, T is the set of technical control replicates, platei is the plate of sample i, and is the raw Olink Ct value of analyte a in sample s. Normalized Ct values ( ′ ) were defined as the z-scores of the plate-adjusted, log transformed Ct values ′ : ′ = ( ) − ( ) Samples with similar normalized Ct profiles ( ′ were clustered using Kmeans++ (https://github.com/tanaylab/tglkmeans) with k=15. Effects of different treatments on Olink cytokine concentrations at different timepoints were calculated by Welch's t-tests with false discovery rate (FDR) correction. Olink protein module scores were calculated by averaging the normalized z-scores of each module analyte. Univariate logistic regression was performed using R package glmnet_4.1. Area under the curve (AUC) values were calculated using R package pROC_1.16.2. Pearson correlation coefficients of Olink clusters were calculated using the averaged values of each Olink analyte per cluster.

CyTOF Data acquisition and analysis
MDIPA stained WB samples were thawed using the SmartTube Prot 1 Thaw/Erythrocyte Lysis protocol. Samples were subsequently barcoded and pooled utilizing the Fluidigm Cell-ID 20-Plex Pd Barcoding Kit and stained with an antibody cocktail against fixation stable markers for more in depth immune profiling. Following sample barcoding and staining, samples were fixed with 2.4% paraformaldehyde in PBS with 0.08% saponin and 125 nM Iridium (Ir) for 30 mins at RT and stored in Cell Staining Buffer until acquisition. Immediately prior to data acquisition, samples were washed once in Cell Acquisition Solution and resuspended at a concentration of 1x10^6 cells/mL for acquisition (including 10% Fluidigm EQ Normalization Beads). The resuspended cells were then acquired on the Helios Mass Cytometer supplemented with a wide bore injector at an event rate below 400 events/second. After data acquisition, samples were de-barcoded using the Astrolabe Diagnostics platform. Cell populations were identified by a combination of an automated approach using the Astrolabe and manual gating as previously described(45).

PBMC preparation for scRNAseq
COVID + patient and COVIDvolunteer PBMC samples were selected based on manual EMR chart review taking into account, and controlling for patient demographics and treatments. Frozen PBMC were thawed at 37°C and resuspended in RPMI-1640 media+ 10% fetal bovine serum (FBS) with 25 U/mL Benzonase before centrifugation at 350 rcf for 5 mins. Cells were resuspended in media and viable cells were counted by AOPI staining in Nexcelom Cellometer Cell Counters. Combinatorial hashes were prepared in wash buffer (PBS + 0.5% bovine serum albumin (BSA)). 500,000 live cells were stained with hashes for 20 mins on ice before 3 washes in wash buffer. Cells were filtered through a 70 μm filter and then a 40 μm filter twice. Filtered cells were counted and loaded with a targeted cell recover of 35,000 cells/lane across 8 lanes of 5' v1.1 NextGEM assay. PBMC samples were processed, hashed, and sequenced in 3 separate batches.

Patient selection for bronchoalveolar lavage
Respiratory samples for research were allocated from BAL obtained from patients ≧ 18 years of age undergoing bronchoscopy with BAL fluid collection for clinical reasons. Patient groups included the following: (1) positive SARS-CoV-2 polymerase chain reaction (PCR) with COVID-19 related acute respiratory failure requiring intubation and mechanical ventilation; and (2) negative SARS-CoV-2 PCR with suspected lung cancer. Patients in Group 2 who had previous positive SARS-CoV-2 PCR or SARS-CoV-2 antibodies were designated as COVID-19 convalescent. Informed consent for bronchoscopy with BAL fluid collection was obtained separately from consent for research. For Group 2, patients were identified by the pulmonologist and physician assistant who performed the bronchoscopy with BAL fluid collection, which was carried out for suspected lung cancer requiring a diagnostic and staging procedure. Negative SARS-CoV-2 PCR was obtained 2 to 5 days prior to bronchoscopy. Patients were selected on the basis of suspected lung cancer not greater than 5 cm in greatest dimension. Patients with known history of, or clinical suspicion for, active infectious or inflammatory lung diseases were excluded.

BAL fluid collection
All respiratory specimens were collected using sterile, flexible, fiberoptic bronchoscopes. Bronchoscopes were flushed prior to the procedure with 5 mL sterile saline, which was collected for research. For Group 1, a single-use bronchoscope was inserted through the endotracheal tube. For Group 2, a reusable bronchoscope was inserted through a laryngeal mask airway or endotracheal tube placed for the procedure. After airway inspection, the bronchoscope was wedged in a distal airway of interest selected by pre-procedure imaging. Sterile saline was instilled in 30 mL aliquots (up to 90 and 210 mL for Groups 1 and 2, respectively) and aspirated. Aspirated BAL fluid was split into parts for clinical use, which included fluid sent for clinical microbiological analysis, and research use, which was transported immediately to the research laboratory on ice and processed as below.

BAL fluid processing
BAL samples were processed within 30 mins of sample collection. Collected BAL was filtered twice through 70 μM filters and centrifuged at 350 rcf for 5 mins at 4°C. BAL supernatant was collected and treated with 0.1% Triton X-100 for 1 hour to inactivate virus before aliquoting into cryovials for storage at -80°C. BAL cells were incubated with Red Blood Cell Lysis buffer (Thermo Fisher Scientific) for 5 mins at RT before washing with PBS+ 0.5% BSA and centrifugation at 350 rcf for 5 mins at 4°C. Viable cells were counted by AOPI staining in Nexcelom Cellometer Cell Counters. Two lanes of 8000 cells from each sample were loaded onto the 10x Chromium Controller for scRNAseq, 5' v1.1 NextGEM assay.

PBMC and BAL scRNAseq data processing
For all scRNAseq datasets, debris and empty droplets were identified with cells that had gene expression (GEX) UMI counts<22. Cells were identified by finding local minimum in the GEX UMI distribution, keeping only cells ≧ the local minimum UMI count. On average, this was 371 UMI counts. For PBMC scRNAseq, hashtag oligo (HTO) UMI counts <22 and excluded from further analysis. For HTO transformation, each feature in the HTO matrix (linear space) was subtracted from its 5 th quantile and divided by its 95 th quantile. Each cell was subsequently divided by its UMI sum. Because hashing reads were not consistently detected, we underwent additional processing to dehash cells.
For initial mapping of cells to samples using only HTOs, a "key" matrix with biological samples as rows and HTO features as columns was created. Each unit in the matrix was populated with a value of "1" if the sample was supposed to be positive for the hashtag and "0" if it was not. For each cell, pairwise distances with a cosine similarity metric were computed from the key matrix to generate a cosine similarity matrix. Samples with the highest cosine similarity could then be assigned to the cell. Initially assigned mappings were called "hto-ini".
Additionally, a signal to noise ratio (SN) was for each cell was calculated by subtracting the highest cosine similarity from the second highest cosine similarity. The SN ratios for each initially assigned sample usually followed a bimodal distribution. The relative minima of the parametric density of SN ratios was used in order identify the local minimum of this distribution. Only cells with SN ≧ the local minimum were kept as the final cells belonging to that sample. We call these final mappings "sample-hto" We used Souporcell to cluster the cells based on the polymorphisms detected from the RNA-Seq alignment(48). We also inputted the genotypes inferred from whole genome DNA-Seq data as a reference in Souporcell in order to map the clusters to respective patients. Next, we leveraged the Souporcell subject mappings along with the initially assigned HTO mappings (hto-ini) in order to deconvolute the patient to the various timepoints. For patients with single timepoints, we assigned the entire Souporcell cluster to the sample. For patients with multiple timepoints, we performed the hashtag de-multiplexing strategy described above to map the cells from the cluster to the patient's timepoints. We called these mapped cells "sample-soc-hto". For final mappings, we intersected the output from both strategies above ("sample-hto" and "sample-soc-hto") in order to get a consensus cell-sample mapping. Only these consensus cells were used for further downstream analysis.

scRNAseq Analysis
Briefly, mRNA reads were tagged with a cell barcode and UMI. These reads were aligned, and count matrices were built. Cell barcodes with at least 500 UMIs were extracted, and cells that were comprised of more than 25% of reads from the mitochondrial genome were filtered from subsequent analysis. The variability in cell counts or UMI counts across samples were not confounding variables in downstream analyses. To account for controlled sampling of each patient sample, each sample was downsampled to 250 cells. The Seurat R package was then used for scaling, batch-aware integration, clustering, dimensionality reduction, and downstream differential gene expression analyses (49-52). To adjust for batch effects, anchors, defined as overlapping shared nearest neighbors, were imputed and used to transform all datasets into a complete shared space. The function SCTransform was used to scale and identify variable genes that constituted principal components for principal component analysis (PCA). The first 15 principal components were used to perform UMAP reduction, once a shared nearestneighbor graph had been generated and clustering had been performed based on the Louvain method for community detection. To identify the cell clusters, cells were down-sampled to 2000 UMIs per cell and variable genes were selected. Gene module analysis was performed by computing a Pearson correlation matrix between genes for each sample using the scDissector R package and grouping highly correlated genes into gene modules by hierarchical clustering(53). Based on this, cell types were manually annotated. Chi-square analysis of gene expression was done across all immune cell clusters, then on specific subsets of lineages, identified subsets of cell states within the myeloid and lymphoid lineages and major cell types.
scRNAseq immune cell cluster frequency correlations, and integrated scRNAseq cell frequencies and Olink proteomics were calculated using the corrplot package (v0.88) and visualized using the pheatmap package (v1.0.12) in R. To identify differentially expressed genes between alveolar macrophages (AM), we computed relative fold change in gene expression between groups and plotted these values with respect to gene expression (total UMI for a gene of interest), represented as a fraction of all UMI per cell. Significant differentially expressed genes were determined using the FindMarkers function from the Seurat R package with default parameters and selecting for genes with an adjusted p-value of less than 0.05.

Lung Autopsy Tissue Section Preparation
Lung autopsy samples were collected within 24 hours of death (average 10.1±6.2 hours) and fixed in 10% neutral-buffered formalin for 24 hours before transfer to 70% Ethanol (EtOH). Samples were then embedded in paraffin and 4 μM tissue sections formalin fixed paraffin embedded (FFPE) sections were cut onto glass slides and baked at 37ºC overnight.
Using an Autostainer (Bond Rx, Leica Biosystems), slides were covered with covertiles (Bond Universal Covertiles, Leica biosystems) and baked for 10 mins at 57ºC. Slides were deparaffinizied in dewax solution and rehydrated in decreasing concentrations of EtOH. Tissue sections were then incubated in antigen retrieval solution (pH 6 or 9) at 95ºC for 20 mins. Tissue sections were incubated in 3% hydrogen peroxide (Bond Polymer Refine Detection Kit DS9800, Leica Biosystems) for 15 mins to block endogenous peroxidases. Next, tissue sections were incubated in serum-free protein block solution (Dako, X0909) for 30 mins to block nonspecific antibody binding. After the first staining cycle, Fab fragments (AffiniPure Fab Fragment Donkey anti-mouse (715-007-003) or anti-rabbit IgG (711-007-003) against that primary antibody species were used to block carry-over staining whenever there was a repeat of same primary antibody species. Primary antibody staining was performed for 1 hour at RT or overnight at 4°C followed by secondary antibody staining. For primary antibody staining, anti-human CD14 (Sigma-Aldrich, AMAB90897) was used at a dilution of 1:1000; anti-human CD68 (Dako, M0814) was used at a dilution of 1:100; anti-human S100A12 (Atlas Antibodies, HPA002881) was used at a dilution of 1:2500; anti-human FABP4 (R&D Systems, AF3150) was used at a dilution of 1:70; anti-human CD66b (BD Pharmingen, 555723) was used at a dilution of 1:600, anti-human Foxp3 (Abcam, ab20034) was used at a dilution of 1:80, anti-human CD3 (Ventana, 90-4341), was used pre-diluted as provided by the supplier; anti-human CD8 (Dako, M7103) was used at a dilution of 1:100; anti-human CD20 (Dako, M0755) was used at a dilution of 1:250. Anti-mouse (Dako K4001), anti-rabbit (Dako K4003), and anti-goat (R&D Systems VC004) secondary antibodies were used pre-diluted as provided by the supplier. Polymer detection system (Bond Polymer Refine Detection Kit, DS9800, Leica Biosystems) was used for horseradish peroxidase signal amplification. Chromogenic revelation was performed using ImmPact AEC (3-amino99ethylcarbazole) substrate (Vector Laboratories, SK4205) for preset incubation times. Slides were counterstained with hematoxylin (Bond Polymer Refine Detection Kit, DS9800, Leica Biosystems).
For manual staining, slides were mounted with a glycerol-based mounting medium (Dako, C0563) and scanned for digital imaging (Hamamatsu NanoZoomer S60 Whole Slide Scanner). The same slides were successively stained, as per MICSSS protocol. Coverslips were removed by placing slides in a rack and immersing in hot tap water at 56ºC until mounting media dissolved. Chemical destaining between stains was performed by immersing slides in gradually diluted EtOH solutions.

MICSSS co-expression analysis
To analyze marker coexpression, a pseudofluorescence composite image of all chromogenic markers was created. The same region of interest (ROI) were selected from images of each marker in QuPath (https://qupath.github.io/). and exported as PNG formatted images without downsampling. Images of different immunostains belonging to the same ROI were transferred to Fiji-ImageJ and co-registered using the TrakEM2 plug-in(54, 55). Color deconvolution was performed using H-AEC vectors for each image to split the RGB images into three 8-bit channels including hematoxylin (blue), AEC (red chromogen color), and residual (green) channel. The best hematoxylin channel was selected as the nuclear channel. AEC channels representing staining of each marker were assigned to different colors by using the lookup tables (LUTs) function of Fiji and hematoxylin channel was assigned to blue color to mimic fluorescent DAPI staining. Next, color inversion was done on all channels and then merged to achieve a multiplexed pseudofluorescent image. We optimized brightness and contrast settings to facilitate visualization for each immunostain channel by comparison with original chromogen images but did not change underlying image pixel values for quantification.

MICSSS quantification
Stained images were scanned at 40x resolution into the .ndpi format and uploaded to Amazon Web Services (AWS) super computer clusters for high-speed analysis on a Python-based Anaconda Jupyter Notebook. Raw red-green-blue (RGB) thumbnail 1.25x resolution images were analyzed using an in-house tissue recognition algorithm that enhanced tissue contours and used optical densities of pooled pixels across the image to determine a tissue mask. The image was then rigidly registered with the corresponding images across all markers for the same tissue within the MICSSS panel. This registration used a third party SimpleElastic package for Python (https://simpleelastix.readthedocs.io/RigidRegistration.html). Following linear registration, the highest resolution image for each marker (40x) was spliced into multiple tiles that spanned approximately 2000 μm in each dimension with 20% surface area overlap in each direction. Each tile was then denoted with an X,Y pooled pixel coordinate value so that the appropriate corresponding tiles across all the markers in the panel for the same tissue would be analyzed together. Each set of tiles was deconvoluted to extract the hemoxylin channel, which remained consistent across all markers. The hematoxylin channel was then registered with an affine registration (which accounts for shear, scale, rotation, and translational dislocation) and a "b-spline" elastic warping to account for any local tissue warping or tissue damage (https://simpleelastix.readthedocs.io/NonRigidRegistration.html) The vector field transformation matrix produced from the high-resolution affine and b-spline registrations was then applied to the raw RGB tile. Registered RGB tiles were analyzed in parallel across the multiple cores of the AWS supercomputer, trimmed to eliminate overlap, and concatenated to produce one final elastically registered RGB image per marker. 100 ROIs of about 500x500 μm were randomly chosen in the image based on where tissue resided. These were chosen from the last tissue mask in the panel to account for any tissue damage or warping. Each of these ROI was processed in parallel across the multiple cores in the AWS supercomputer. Next, we used the Stardist package for Python (https://github.com/stardist/stardist) that was trained with hematoxylin and eosin staining, to segment each ROI, and to determine the centroids and morphological properties of each determined cell. We previously optimized the sensitivity of this algorithm with these tissues and therefore used an overall sensitivity value of 0.1 for the algorithm. This algorithm provides information for the nucleus of each cell in the ROI, which we then artificially expanded by 5 μm to simulate the cell's cytoplasm. Overlapping cytoplasm from adjacent cells in dense regions were condensed to prevent overcounting surface area.
Each cell per ROI was then analyzed for marker positivity. ROI were deconvoluted for its AEC detection channel and each cell was translated to a median AEC value from pixels that resided in its nucleus and artificially expanded cytoplasm. If the median AEC value for the cell was above the threshold value for positivity deemed for that marker, the cell was considered positive. Threshold values per marker were determined with a pathologist. The percentage of positive cells was determined by the number of cells above the threshold for that marker over the total number of cells in that ROI. Co-expression analysis was performed by determining if a cell was positive for the multiple markers of interest over the total number of cells in that ROI. We plotted each value per ROI.  TRAC  CD8A  CD8B  CD4  CD40LG  FOXP3  CCR7  IL7R  LDHB  LT B  PDCD1  CTLA4  TIGIT  LAG3  TRDC  KLRC1  XCL2  LITAF  ID2  AREG  CD7  TYROBP  FCER1G  KLRF1  FCGR3A  FGFBP2  KLRB1  ALOX5AP  KLRD1  NKG77  GZMK  GZMB  GNLY  CTSW  GZMA  CCL5  GZMH  GZMM  CCL4  HLA-DRB1  FA S  CD27  CD28  TBX21  CD38  IL32  CX3CR1  LEF1  TCF7  CCL3 CD79B  TFRC  SDC1  CD38  CD27  IGHG1  IGHG3  IGHG4  IGHD  IGHM  IGHA1  IGHA2  CALR  HSP90B1  PPIB  SEC11C  JCHAIN  PRDM1  XBP1  MZB1  HLA-DQB1  HLA-DRA  CD74  TXNIP  FCER2  FCMR  SELL  BANK1  ITGAX  CD69  DUSP2  TNFRSF13B  PTPRC  CXCR5  NT5E  COVIDsamples were obtained from healthy volunteers (A to G). Statistical significance in (A, B, D, E, and G) were determined by Kruskal-Wallis followed by multiple comparisons test with false discovery rate correction; ns, not significant; *q<0.05; **q<0.01; ***q<0.001; ****q<0.0001.